!----------------------------------------------------------
!----------------------------------------------------------
subroutine qmga(step)
  use rbf_type
  use variable, only: nodes,a,beta,nn,xm,wall,npatch
  use patch_module, only:surface_nano,inside_nano,pm
  implicit none
  integer (i4b) ,intent(in) :: step
  integer (i4b) :: i,ii,part
  real(dp),dimension(3) :: e,e_ref,rotvec,director
  real(dp) :: qw,qx,qy,qz,pdot,norm_a,norm_b,angle,jj,kk,mm,norm_vec,scalar
  real(dp), dimension(5) ::   q_tmp
   integer (i4b) ::  iout=50
!*********************************************************************
  Character(len=10) :: namefile
  Character(len=20), Parameter :: sfile_1= 'traj.'
  
  write(namefile,'(I10)')step
  open(iout,FILE=trim(adjustr(sfile_1))//trim(adjustl(namefile)),STATUS='replace')

!**************************************************************************************       
  
  write(iout,100)'ITEM: TIMESTEP'
  write(iout,102)iout
  write(iout,100)'ITEM: NUMBER OF ATOMS'
  write(iout,102)nodes
  write(iout,100)'ITEM: BOX BOUNDS'
  write(iout,101)-beta*0.5,beta*0.5
  write(iout,101)-beta*0.5,beta*0.5
  write(iout,101)-beta*0.5,beta*0.5
  write(iout,100)'ITEM: ATOMS id type x y z c_qdat[1] c_qdat[2] c_qdat[3] c_qdat[4] '
 
  e_ref(1)=0.d0
  e_ref(2)=0.d0
  e_ref(3)=1.d0


  do i = 1,nodes
     ii=i*3
     call inv_transformation( a(:,i), q_tmp(:))
     call eigenval( q_tmp(:), scalar,director(:))

     e=director(:) 
     norm_a = sqrt(e(1)**2 + e(2)**2 + e(3)**2)
! Normalize the vectors

     e(1)=e(1)/norm_a
     e(2)=e(2)/norm_a
     e(3)=e(3)/norm_a

     if (e_ref(1) == e(1).and.e_ref(2)== e(2).and.(e(3)==e_ref(3).or.e(3)==-e_ref(3)))then
        rotvec(1) = 1.d0
        rotvec(2) = 0.d0
        rotvec(3) = 0.d0
     else 
        rotvec(1) =   e_ref(2)*e(3)- e_ref(3)*e(2)
        rotvec(2) = -(e_ref(1)*e(3)- e_ref(3)*e(1))
        rotvec(3) =   e_ref(1)*e(2)- e_ref(2)*e(1)
     endif

     norm_vec = sqrt(rotvec(1)**2 + rotvec(2)**2 + rotvec(3)**2)

     rotvec(1) = rotvec(1)/norm_vec
     rotvec(2) = rotvec(2)/norm_vec
     rotvec(3) = rotvec(3)/norm_vec

     angle = acos(e(1)*e_ref(1) + e(2)*e_ref(2) + e(3)*e_ref(3))

     qw = cos(angle/2.d0)
     qx = rotvec(1) * sin(angle/2.d0)
     qy = rotvec(2) * sin(angle/2.d0)
     qz = rotvec(3) * sin(angle/2.d0)
     

        if(i.le.nn)then
           if(surface_nano(i))then
              write(iout,458)i,3,xm(ii - 2:ii),qw,qx,qy,qz,' '
           else if (inside_nano(i))then
              write(iout,458)i,4,xm(ii - 2:ii),qw,qx,qy,qz,' '
           else
              write(iout,458)i,1,xm(ii - 2:ii),qw,qx,qy,qz,' '
           endif
        else
           if(surface_nano(i))then
              write(iout,458)i,3,xm(ii - 2:ii),qw,qx,qy,qz,' '
           else if (inside_nano(i))then
              write(iout,458)i,4,xm(ii - 2:ii),qw,qx,qy,qz,' '
           else
              write(iout,458)i,2,xm(ii - 2:ii),qw,qx,qy,qz,' '
           endif
        endif

  end do
  if(npatch.ne.0)then

     do part=1,npatch
        write(iout,458)nodes+part,5,pm(3*part - 2:3*part),1.0,0.0,0.0,0.0,' '
     enddo
  endif

  close(iout)
458 format(I10,1x,I3,1x,7(f14.7),A)        
  
100 format(A)
101 format(2f7.2)
102 format (I8)
end subroutine qmga

!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------

subroutine paraview_eigen_new(step)
  use rbf_type
  use variable, only: nodes,a,Nx,Ny,Nz,beta,xm,inside,surface,nee_cube,ddx,ddy,ddz,xmax,ymax,zmax,ry,rz
  use patch_module, only:surface_nano,inside_nano,pm
  implicit none
  integer (i4b) ,intent(in) :: step
  
  integer (i4b) :: i,ii,iout=50,j,k,iii,count_mol
  real(dp), dimension(5) ::   q_tmp,tensor
  real(dp) :: r2,lbox,cx,scalar,x,y,z,eta,factor,subtract
  real(dp),dimension(3) :: dx,director,xs,dir_vector
!*********************************************************************
  Character(len=10) :: namefile
  Character(len=20), Parameter :: sfile_3= 'isosurface.',end='.vtk'  
  
  write(namefile,'(I10)')step
  open(iout,FILE=trim(adjustr(sfile_3))//trim(adjustl(namefile))//trim(adjustl(end)),STATUS='replace')
  
  xs(1) = 0.5d0*xmax
  xs(2) = 0.5d0*ymax
  xs(3) = 0.5d0*zmax

!**************************************************************************************

  write(iout,100)'# vtk DataFile Version 3.0'
  write(iout,100)'vtk output'
  write(iout,100)'ASCII'
  write(iout,100)'DATASET STRUCTURED_GRID'
  write(iout,110)'DIMENSIONS',Nx,Ny,Nz
  write(iout,120)'POINTS', Nx*Ny*Nz ,'  float'
  
  do iii=1,nz
     do j=1,ny
        do k=1,nx 
           
           x = ddx*(k-1)
           y = ddy*(j-1)
           z = ddz*(iii-1)
           write(iout,151)  x - xs(1),y - xs(2),z - xs(3),'  '
           
        enddo
     enddo
  enddo


  write(iout,*)
  write(iout,120)'POINT_DATA',nx*ny*nz
  write(iout,100)'SCALARS scalars float'
  write(iout,100)'LOOKUP_TABLE default'

  do k=1, nz
     do j=1, ny
        do iii=1, nx

       
              i=nee_cube(iii,j,k)
              call inv_transformation( a(:,i), q_tmp(:))
              call eigenval       ( q_tmp(:), scalar, director(:))
              write(iout,152 )  scalar*1.5d0,'  '
     
              
        enddo
     enddo
  enddo



count_mol = 0


    do k=1, nz
     do j=1, ny
        do iii=1, nx
           
              i=nee_cube(iii,j,k)
              if(.not.surface_nano(i).and..not.inside_nano(i))then
                 count_mol = count_mol + 1
                 call inv_transformation( a(:,i), q_tmp(:))
                 call autovalores_biax(q_tmp(:), scalar, eta,director(:))
                 

                 tensor (1)= tensor(1)+ director (1)*director(1)
                 tensor (2)= tensor(2)+ director (1)*director(2)
                 tensor (3)= tensor(3)+ director (1)*director(3)
                 tensor (4)= tensor(4)+ director (2)*director(2)
                 tensor (5)= tensor(5)+ director (2)*director(3)
                
             endif
           
        enddo
     enddo
  enddo

  factor=3.d0/(2.d0*count_mol)
  subtract = 1.d0/2.d0

  tensor (1)= factor*tensor(1)-subtract
  tensor (2)= factor*tensor(2)
  tensor (3)= factor*tensor(3)
  tensor (4)= factor*tensor(4)-subtract
  tensor (5)= factor*tensor(5)

  call autovalores_biax(tensor(:), scalar, eta,dir_vector(:))


  write(iout,*)
  write(iout,100)'SCALARS dir_vector float'
  write(iout,100)'LOOKUP_TABLE default'

   do k=1, nz
     do j=1, ny
        do iii=1, nx
         
              i=nee_cube(iii,j,k)
       !       if(.not.surface_nano(i).and..not.inside_nano(i))then
                 
                 call inv_transformation( a(:,i), q_tmp(:))
                 call autovalores_biax(q_tmp(:), scalar, eta,director(:))
             
                 
                 
                 write(iout,152 )abs(dot_product(director,dir_vector))
      
       !       endif
        enddo
     enddo
  enddo


  write(iout,*)
  write(iout,100)'VECTORS directors float'

  do k=1, nz
     do j=1, ny
        do iii=1, nx
      
              i=nee_cube(iii,j,k)
              call inv_transformation( a(:,i), q_tmp(:))
              call eigenval       ( q_tmp(:), scalar, director(:))

                 if(director(3).lt.0.d0)then
                    director= -1.d0*director
                 endif
              write(iout,153 )  director,'  '
      
              
        enddo
     enddo
  enddo
  

  close (iout)
!*****************************************************************************************
!*****************************************************************************************
100 format('',A)
120 format('',A,I12,A)
110 format('',A, I6, 1X, I6, 1X, I6)
151 format (3f12.6 ,A)
152 format (f12.6 ,A)
153 format (3f12.6 ,A)
154 format (f12.6 )
155 format (3(f8.4,2x) )
end subroutine paraview_eigen_new

